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ABSTRACT 

Recent observations of strikingly well-defined spirals in the circumstellar envelopes 
of asymptotic giant branch (AGB) stars point to the existence of binary companions in 
these objects. In the case of planet or brown dwarf mass companions, we investigate the 
observational properties of the spiral-onion shell wakes due to the gravitational interac- 
tion of these companions with the outflowing circumstellar matter. Three dimensional 
hydrodynamical simulations at high resolution show that the substellar mass objects 
produce detectable signatures, corresponding to density contrasts (10-200 %) and arm 
separations (10-400 AU) at 100 AU distance, for the wake induced by a Jupiter to brown 
dwarf mass object orbiting a solar mass AGB star. In particular, the arm pattern prop- 
agates in the radial direction with a speed depending on the local wind speed and sound 
speed, implying possible variations in the arm separation in the wind acceleration region 
and/or in a slow wind with significant temperature variation. The pattern propagation 
speeds of the inner and outer boundaries differ by twice the sound speed, leading to 
the overlap of high density boundaries in slow winds and producing a subpattern of 
the spiral arm feature. Vertically, the wake forms concentric arcs with angular sizes 
anticorrelated to the wind Mach number. We provide an empirical formula for the peak 
density enhancement as a function of the mass, orbital distance, and velocity of the ob- 
ject as well as the wind and local sound speed. In typical condition of AGB envelopes, 
the arm-interarm density contrast can be greater than 30 % of the background density 
within a distance of ~ 10 (Mp/Mj) AU for the object mass Mp in units of Jupiter mass 
Mj. These results suggest that such features may probe unseen substellar mass ob- 
jects embedded in the winds of AGB stars and may be useful in planning future high 
sensitivity /resolution observations with ALMA. 
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INTRODUCTION 



In high density environments, substehar companions leave their orbital imprints through the 
gravitational interaction with the surrounding circumstellar mediur n. An example is the ob servation 
of molecular spiral arms in the protoplanetary disk of A B Aurigae (IFukagawa et al.ll2004l '). which is 
considered as evidence of a giant planet within the disk (ILin et al.ll2006l). It is now well established 
that a spiral wave sets up in the differentially rotating system (e.g., iGoldreich &: Tremainelll979l : 
Massetll2008l . and references therein) during the formation and growth of a planet. This protoplan- 
etary phase lasts for a few million years and is followed by a main sequence phase, in which the 
circumstellar disk dissipates and the planet-disk interaction ceases. 

Planets or brown dwarfs orbiting low and intermediate mass stars {M^, < 6 — 8Mq) may have 
another opportunity to interact with their immediate surroundings before the stars terminate their 
evolution as white dwarf remnants. Stars that have evolved off the main sequence, especially on the 
asymptotic giant branch (AGB), expel a large fraction of their mass at rates ~ 10^^ — 10"^ Mq yr~^, 
forming extensive (10^^ — 10^^ cm) cool (10^ — 10^ K) low velocity (3— 30kms~^) outflowing envelopes 



(jHabing Sz OlofssonI l2003l : iFong et al.l l2006l . and references therein) . Such intense mass loss rates 
imply an envelope density of 10~^^ 



10 gem ^ at 100 AU dista nce from the stellar center, whic h 



is slig htly lower than the typical density of protoplanetary disks (jllavashil Il98ll : lAikawa &: Herbst 
19991 ). Similar to the protoplanetary case, substellar mass objects that are embedded in the dense 
circumstellar envelopes of evolved giant stars will interact with the environment, playing a role in 
shaping the wind structure and also affecting the orbital evolution of the system. 

Recently, the circ umstellar envelopes of a few AGB stars h ave been found to possess spi- 
ral morphologies (e.g., iMauron &: Huggind l2006l for AFGL 3068; iDinh-V.-Trung &: Liml l2009l for 
CIT 6). Such structures cannot be expl ained by the spherically sy mmetric wind mod els for 
isolated stars with either continued (e.g.. iParkerl 119581: iGilmanI Il972l ) or pulsating (e.g., IWood 



1979 



Willson fc Hilllll979l : IWinters et al. 



200C 



Simis et al.ll200ll ) mass ejection (for a review, see 



Lafon fc Berruveij|l99ll ). although the pulsation models are able to explain the conc entric shells 



and arcs of dust detected in a number of 



1998 



Mauron &: Hugging 



199S 



Kwok et al. 



30st-AGB or proto-plane tary phases (e.g., ISahai et al 



2001 



Hrivnak et al. 



2001 



However, direct comparison 
with the detailed observational aspects remains as they likely involve the nonlinear process of dust 
formation in the envelopes. On the other hand, the spirals and greatly distended shells/arcs are 
more naturally produced by dynai nical effects arising from a binary interaction in the expanding 



wind of the mass losing AGB star (jSokeiill994l : iMastrodemos &: Morris 



199S 



Hel l2007l : lEdgar et al 



20081 ). The binary- induced spiral model is supported by t he detection of two point sources within 
the amazingly well-defined spiral envelope of AFGL 3068 ([Morris et al.l 120061 ). 



Interestingly, albeit initiated from a differe nt viewpoint, a very similar spiral-onion shell struc- 
ture in three dimensional space was described by lKim fc KimI (j2007l ). This study may help facilitate 
an understanding of asymmetric morphologies in circumstellar envelopes for cases of companions 
with low masses, which are not sufficiently massive to displace the parent star due to orbital motion. 
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Their linear perturbation analysis f or gravi tationally induced density wakes of circularly orbiting 
objects was motivated by lOstrikerl 's (|l999l ) purely analytical formali sm of the densit y wakes for 
linear orbit objects. Follow up st udies were carried out for bi naries (jKim et al.ll2008l ). for heavy 



objec ts showing nonl in ear aspects (IKim &: Kimll2009l : lKimll201[ll ). and for objects in accelerated mo- 



tion (jNamounilboiol ). kim k KinJ ^20od ) noticed that the local features of the nonlinear density 
wake are conceptually the same as those of a Bondi-Hoyle-Lyttleton accretion column (see lEdgar 
2004I . and references therein), and are insensitive to the surface condition of the object. 



We note that the above theoretical works are based on an initially static uniform medium, 
which differs from the expanding AGB envelopes. However, their results deserve to be examined 
in the extreme limit of a slow wind. The focus of these works was on the orbital evolution of the 
perturbing obje ct due to d ynamical frict ion, leaving the details of the wake itself in abeyance. In a 
separate paper ( Kim|[2011 . hereafter Kll ) the properties of the wake induced by a circularly orbiting 
object at distance with the orbital Mach number Aip (>1) in a static uniform background was 
investigated. Here, we compare the results for the outflowing envelopes to the case of a static 
background. The resulting wake in the latter is characterized by a single-armed Archimedes spiral 
in the orbital plane satisfying r/vp = ipM~^ + 1, where Lp is the angular distance along the spiral 
pattern, and by arcs in meridional planes having the angular size 

of 2tan~i(7W2-l)V2. its density 

contrast is determined by 

TB 1 



{Ml 



(1) 



1)1/2' 

and the background density modification due to memory effect, where is the Bondi accretion 
radius of the perturbing object. In this paper, we investigate the effects of an outflowing back- 
ground on the wake features both analytically and numerically in order to provide an interpretative 
framework of future observations for indirectly probing substellar mass objects in AGB envelopes. 
In ^ we describe the setup of the numerical simulation by specifying the stellar wind model as a 
background and the gravitational perturbation due to an object orbiting about the central star. In 
^we present the numerical results focusing on the shape and morphology of the wake as well as 
its density contrast to the background matter. Finally in ^ we discuss our findings with a view 
toward future observations. 



2. SIMULATION SETUP 

2.1. Background Wind Condition 

For a wind characterized by a velocity of V^, the density distribution follows from the 
steady hydrodynamic conditions 

V • (p^Vw) = 0, (2) 
and ^ 

Vw • VVw = - ^ Vp^ - V$f , (3) 

Pw 
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for a background gravitational potential ^'f^ . We do not explicitly consider the effects of dust, 
which is believed to play an essential role in driving the wind of cool giant stars, but we use an 
effective potential term by introducing an effective stellar mass Mf = M*(l - /) conceptually 
to include an additional force due to stellar radiation pressure on the wind, //r^. Here, Cg represents 
the sound speed of the background gas flow, = ^p/p with the adiabatic index 7 in general. In 
the numerical simulations, we employ an isothermal gas approximation (7 = 1). 

In a spherical description, the wind material has a density distribution = /9to(r) stratified 
in the radial direction with a corresponding velocity field = V^(r)f affected by the pressure 
gradient. Given a continuous stellar mass loss rate M*, mass conservation (eq. [2]) for a radial 
wind leads to Pw{f) = M^/Anr'^Vwir). A constant velocity, as is observed at large radii, results in 
the wind attaining a density distribution. The details of the velocity structure are regulated 
by the driving mechanism of the wind. For instance, a purely gaseous wind without any additional 
force in the momentum equation (eq. [3]) approaches a near-constant velocity regime only at large 
distances whe re both the stellar gravity and the gas pressure gradient have significantly decreased. 
Parked (jl958l ) derived the isothermal wind solution as 



Ml - Ml{r,) - 2ln 



41n( -) +4(^ 



(4) 



where Ad^j = V^/cg represents the wind Mach number and the critical radius = GM^/lc^ is the 
so-called sonic radius in a transonic wind solution satisfying Mwifc) = 1- A hydro dynamical wind 
from the star follows the supersonic (A^^(rc) > 1), transonic {^Awirc) = 1), or subsonic (Al,,,(r^) < 



1) branch of iParkeij's wind solutions (Fig. [T^-b; see also Figure 3.1 of lLamers &: Cassinellilll999l ). 
Any branch of iPar kerf s wind solutions shows significant variation of the outflowing velocity inside 
the critical radius on the scale of a few tens to thousands AU for typical red giant stars (see below 
for accelerating winds). 

The additional outward force due to the stellar radiation pressure on dust , which couples with 
the gas, is commonly rega rded as the wind driving mechanism of AGB stars (jLamers Cassinelli 



199S 



Winters et al.l 120001 ) accelerating the wind faster so that a near-constant velocity is reached 
at smaller radii. Strong acceleration (/ = 1; = 0) modifles the isothermal wind to satisfy (see 

Fig. m 



Ml-Mlin)-2ln 



41n ( — 



(5) 



whose quantities are determined by the conditions at the dust condensation radius r^. In reality, 
this is probably not a distinct boundary as the wind accelerates gradually and / would not reach 
unity. Our approximations oversimplify the wind driving mechanism, but'they are adopted for the 
background wind condition in order to focus on the larger scale structure. 

Although the background wind condition depends on the wind driving mechanism, the response 
of gas to the gravitational perturbation is insensitive to the global wind structure, but instead is 
determined by the local wind quantities. This point will be checked by two sets of simulations 
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characterized by the transonic and supersonic branches of iParkeij 's wind solutions, as weh as one 
set of simulations with / = 1 mimicking strong radiation pressure. Hence, it will be confirmed 
that the wake properties are all determined by the wind density Pwir), its expansion speed T4,(r), 
and the sound speed Cs(r) on the spot in addition to the orbital and accretion properties of the 
perturbing substellar mass object. 



2.2. Model Parameters 



In the background wind characterized by its outflowing velocity and corresponding density 
Pw, we introduce a perturbation (5V for velocity and 5p for density due to the gravitational potential 
of a substellar companion of mass Mp orbiting about the central star at a distance with the 
velocity Vp. The response of the gaseous medium is governed by the ideal hydrodynamic equations 



dp ^ 



(pV) = 0, 



and 



av 

'dt 



+ v- vv 



p 



-Vp - v^> 



cff 



(6) 



(7) 



for the density p = pw + Sp and the velo c ity fi elds V = + 5V. These equations are inte- 
grated by the FLASH code of iFrvxell et al.l ([20001) in three dimensional Cartesian coordinates. We 
utilize the adaptive rn esh refinement (AMR) implementation based on the PARAMESH package 
(jMacNeice et al.lll999l ). in order to minimize the deviations of the wind from a spherical shape 
by assigning sufficient number of grid cells within from the center of the mass losing star lo- 
cated at the origin of the Cartesian coordinates. However, to avoid potential artifacts due to the 
rapidly propagating spiral shocks, the refinement is set only at the beginning of each run such 
that the nonuniform subgrids are stationary through out the simulation. We use the maximum 
level of refinement up to eight with the mother grid number of 64x64x32. For the calculation of 
total gravitational field, the POINTMASS implementation is used with small modifications (1) to 
include two sourc es of gravity, (2) to introd uce the gravitational softening radii of Plummer-type 
objects (see, e.g.. lBinnev &: Tremaind 120081 ) . and (3) to allow the motion of the objects. 



As the perturber, we consider a substellar mass object with the mass Mp greater than 10 
Mj where Mj denotes the mass of Jupiter. The orbit of the perturber is assumed to be circular, 
ignoring evolutionary processes such as stellar mass loss, mass transport to the substellar object, 
tidal interaction between th e star and the object, and frictional and gravitational drag forces (see. 



e.g., iVillaver &: Livid l2009l . and references therein). The range of orbital radius rp considered in 
our simulations is 10-150 AU. We also confirm that the effect of the perturber size, defined by its 
gravitational softening radius r^, is minor when its size is sufficiently small, exponentially decreasing 
with decreasing (0.1-10 AU). The key parameters are the velocity of the object, Vp, and of the 
wind in the outflowing motion Vw, relative to the sound speed Cg. In the numerical simulations, we 
treat only a pseudo-isothermal gas (7 = 1.00001) with the constant sound speed Cg ranging between 
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1-10 kms~ . Given Cg, the object speed is arbitrarily selected in the range of Mp = 0.5-10, to 
extend parameter space. The wind speeds correspond to a Mach number Aiy^ < 11. We perform 50 
models in total, and Table [J lists the model parameters of our selected sets of simulations and the 
respective resolutions. The Bondi accretion radius of the perturber, rs = GMp/c^, is also listed 
for reference. 

The simulation sets are categorized into accelerated (A), supersonic (S), and transonic (T) 
background wind models. The S models are further divided according to the wind Mach number 
into the fast {Mw = 10) SF and slow (TW^ = 5) SS models. For the A models that treat the wind 
accelerated by an arbitrary force, we focus on highly supersonic winds, in which the velocity 
variation is significantly less than that of transonic or marginally supersonic winds, facilitating 
analysis of the effects of the other parameters explicitly related to the orbiting object. For the S 
and A models, we start with a constant outflow velocity and a density distribution in the entire 
domain, and reset the quantities inside every timestep to avoid the formal infinite values at the 
center (eqs. [1] and [5]) if an inner boundary is not set. Eventually the flow follows the analytic 
solution, either equation dl]) or ([5]), determined by the conditions at the reset radius r^, (Fig. [T]). 
The background density py^ and velocity field for the S and A models are obtained from the 
separate runs in which the gravitational influence of the perturbing object is suppressed. 



3. RESULTS 

Figure [2] shows the density enhancemen t fact or, a = 6p/pw, for a fast wind model (Al; see 



Tabled]). Compared to the wakes analyzed in lKllI excluding the background wind, the wake in the 



outflowing environment is radially more extended by it 1. It is also concentrated to a greater 
extent toward the orbital plane showing a correlation with both wind and object Mach numbers, 
and A4p. In this section, we elucidate the details of the shape and morphology of the wake 
pattern in the outflowing circumstance, and determine the arm-interarm density contrast and its 
dependence on the wind and object properties. 



3.1. Shape of Spiral Pattern 

The shape of the wake is highly affected by the background flow when formed by the gravi- 
tational interaction of the orbiting object with the background medium. In particular, the wake 
in the wind of an AGB star forms a loosely wound spiral in the orbital plane in contrast to the 
case of a static medium. The shape of the spiral is determined by the ratio between the pattern 



propagation speed 14,rm and the orbital speed of the perturbing object Vp (see lKll 



d{r/rp) _ Varm 
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which reproduces T4rm = Cs in lKim fc KimI (j2007l ) when the wind is absent. In Figure [3l we compare 
the spirals induced by the gravitational potential of the orbiting object in different background 
velocity fields. The change of the spiral shapes due to the different wind velocities is highlighted 
by the dashed line showing the opening angle of the spiral in a static medium. This line matches 
the spiral shape in the vicinity of the object in Figure [3^, while the spirals open further with 
higher wind speeds as displayed in Figure [3b-c. This suggests that, in the presence of wind, the 
propagation speed of the spiral pattern in the radial direction, T4mn is faster than the sonic speed 
of the medium. 

For a linear perturbation analysis, we assume that the wind velocity is radial and constant 
in the entire region (see i j2.1l for the background wind conditions applied in numerical simulations). 
Combining equations dS])-© in the linearized versions for the perturbed velocity (5V under the 
background conditions satisfying equations ©-([S]), we obtain 



D^SY 



V 



-f I • (5V 

r 



Dt ^ 



(9) 



where D/Dt = d/dt + \^ ■ V is the time derivative in the frame moving with the local background 
wind. Due to the external gravity on the right hand side of the equation, the perturbed gas 
forms a single - armed spiral-like structure, similar to one in the static background analyzed in 
detail by 



Kll 



/ (Kr)2-(fcr)2 


i{kr + i) {A + cot 9) 


-{kr + i)M 


\ 


(ir \ 




iAkr 


{Krf - L2 


iM{A - cot 6 


) 


= 


{') 


-Mkr 


iM{A + cot 6) 


{Krf - M2 


/ 







Using spherical harmonics for the perturbation, in general, in the form of S\ = 
i^(r)Pj™'(cos 9) exp{ikr + imip — icot), where ^ is a slowly varying function of radius, we can write 



(10) 



where K, A, M, and are defined as K = {oo — T4)/c)/cs, A = In P™ (cos 6*), M = m/ sin 9, and 
= l{l -|- 1) -|- (1 — m^)/sin^0, respectively. In the WKB approximation {kr ^ 1 and Kr ^ 1), 
the dispersion relation yields to = {Vw ± Cg) k, indicating that the propagation speed of the pattern 
is lb Cs at large distances, where and Cs are local quantities. 

In order to check T4rm in the numerical simulations, a density map is displayed in Figure [4^ 
for the pattern in polar coordinates ((/?, r) in which the slope of the high density arm boundaries 
represents the velocity ratio Varm/Vp- Numerical differentiation along the high density boundary 
indeed gives the propagation speed of V^, -|- Cs for the outer boundary and — Cs for the inner 
boundary (Fig. Hb), depending on the local wind properties and not on the background wind type. 
This result implies that in the warm central region where the thermal sound speed is relatively 
large, the shape of spiral must take account of the sound speed as well as the expansion speed 
of the envelope. The spacing of the arm is not necessarily constant, but varies especially in the 
wind acceleration zone and in slow wind envelopes. In Figure the outer boundary reaches 
r/rp — 1 = 2.0 and 4.3 after the first and second turns (see also vertical dotted lines in Fig. [JJ^) 
thus the arm spacing increases by 15% (from Ar/r^ = 2.0 to 2.3). For the inner boundary, the 
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locations of r/r^ -1 = 0, 0.5, 1.2, 2.0, 3.0, 4.2, and 5.4 indicate Ar/rp = 0.5, 0.7, 0.8, 1.0, 1.1, and 
1.2 for the respective turns. That is, the arm spacing changes significantly at the first a few turns, 
where the wind velocity varies considerably in this model, and the increase rate reaches to 140% 
at the sixth turn relative to the arm spacing at the first turn. 

Of particular importance is the overlap of the inner and outer arm boundaries seen in Figure Hti 
at Tover/fp = 1.7, 3.0, and 5.1. In outfiowing envelopes the outer and inner boundaries propagate 
at different s peed s (by 2cs), causing overlaps, unlike the unique propagation speed Cg in a static 



background (jKllI ). Thus the detailed structure of the spiral in the outflowing envelope can be 
more complicated than in a static background. The difference in propagation speeds results in 
the overlap of boundaries at resonant positions rover /'"p given by the integration of equation ^ as 
iVy, + Cg) /Vp X ~ iy^j — Cs) /Vp X + 27r X integer), assuming that the pattern propagation speed 
does not change significantly. Hence, the overlap occurs every Arover derived as 

Tp Mp 

consistent with the intervals between overlays in Figured^. 

We should also notice the presence of spiral structures in the cases of subsonic objects (SFx 
and SSx models). In a static background medium, an object with a subsonic orbital speed {Vp < 
Cs) induces the wa ke to be in a smooth distribution rather than confined within a spiral arm 



(jKim &: KimI 120071 ). However, in an outfiowing environment, the background wind sweeps and 
compresses the induced wake material, resulting in a spiral arm shape, as in the cases for supersonic 
motion. 



3.2. Density Enhancement 

Figure [5^ exhibits the density enhancement factor for model Al. In order to trace the full 
structure of the spiral, we plot the profiles along eight different directions in the orbital plane with 
different colored lines. One noticeable difference from the corresponding profiles in a static back- 
ground is the level of the zero baseline of the perturbed density. Without a wind, the wake always 
has minimum density with the central value of a = rs/vp forming a hydrodynamic equilibrium (see 



Kll 



for details). However, in the outfiowing medium, this global increase of background density 
is suppressed since the material is swept out of the central part by the wind before its pressure 
balances the gravitational potential at the central part. This wind sweeping effect also precludes 
the wake from developing the nonlinear features described by iKimI ( 20ld ). In Figure [5^1., the peak 
values of the density enhancement decrease with distance as outlined by a solid line (see below). 
This is approximated by a dashed line, a = (r^jr — rp\~^)^^'^, with the power law index of 0.5 
except for the steeper profile in the innermost part. 

In order to explore the dependence of the density enhancement a on the properties of object 
and envelope, we compare in Figure [5] the profiles of a along distances in the orbital plane for the 
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A models. Compared to the fiducial model (Al), we investigate the effects of the stellar mass loss 
rate M^, (A2), the size of the star defined by the numerical reset radius (A3), the gravitational 
softening radius of the perturbing object (A4), the mass of the object Mp as well as the sound 
speed of the background gaseous medium Cg (A5), and the orbital radius of the object rp (A6). By 
comparing model Al with models A2 and A3, it is found that the perturbed density a = 5p/pw 
does not depend on the mass loss rate and the radius r* at which the mass loss occurs, while 
the overall density structure Pw{f) = M*/47rr^Vu,(r) is regulated by the mass loss rate Af* and 
the background wind velocity The independence of density enhancement on the mass loss 

rate is a natural consequence for an isothermal gas. Reducing the gravitational softening radius 
r<j of the perturbing object by an order of magnitude from 0.1 (Al) to 0.01 (A4) sharpens 
the density profiles to some extent but insignificantly changes the overall features. The identity of 
density profiles in the models A4 and A5 indicates that the density enhancement a depends merely 
on the Bondi accretion radius tb oc Mp/c^, and not on the object mass Mp or the sound speed Cg 
individually; it is also checked that a is proportional to Mp and c^^ using models A5-1 and A5-2, 
respectively. For model A6, the perturbing object is located at a distance twice as large from the 
central star than the orbital distance in the other A models, resulting in a similar value of a at 
large distances r ^ rp although the width and interval of the spiral have increased with the larger 
orbital radius rp. 

In Figure [H we further examine the physical parameters that affect the wake density using the 
S models, focusing on the velocities related to the object motion (Vp) and the wind expansion (Vw) 
in units of background sound speed Cg. Comparison between SF models (Fig. [6H-c) and between 
SS models (Fig. [6li-f) separately reveals a definite anticorrelation of the pattern spacing with the 
orbital Mach number Mp given a nearly constant wind Mach number A4w (9.5-10 in SF; 3.9-5.2 
in SS). On the contrary, one can also find the decrease of the pattern spacing with a smaller wind 
Mach number M^i by comparing the profiles in (a) and (d), (b) and (e), or (c) and (f). The reduced 
spacing of the spiral arm pattern facilitates overlaps between the inner and outer arm boundaries 
characterized by different propagation speeds with respect to each other. As a consequence of 
the overlaps, the wake develops subpatterns showing additional increases in peak densities at the 
overlap positions, rover j with gradual decreases between the overlap positions. A rough estimate 
for the overlap intervals, Arover/fp from equation (jlip . with the assumption of constant A4w is 
consistent with the intervals between the highest density peaks in Figured > 50 for (a) and (b), 
~ 30 for (c), 20 — 30 for (d), 10 — 15 for (e), and 5 — 8 for (f). The peak densities in Figure [6] 
(excluding the extra peaks due to the subpattern) tend to increase with lower orbital Mach number 
and higher wind Mach number at r ^ rs, albeit the dependence on the wind Mach number is not 
as clear. 

We generalize all the parameter dependence described above, using 50 models of A-, S-, and 
T-type background winds in the range of A4w < 11 with Cg = 1-10 kms~^ for p erturbing objects 



of Aip = 0.5-10 and rs/rp = 0.05-1. By modifying equation ([T]), based on lKllP s analysis for the 



simpler case of a static uniform background, we obtain an empirical formula for the minimum value 
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of the density peaks (black solid lines in Fig. [5] and [6]) described as follows: 

which yields ai (eq. [1]) when Aiy^ = 0. We note that this equation is suitable for supersonic winds 
since the winds in our simulations are subsonic only in the innermost region in T-type models, 
which does not sufficiently contribute to the overall profile included in our analysis. A subsonic 
wind likely leads to a more complicated wake description because (1) it is not able to fully sweep 
off the delayed perturbations from the object in the past orbits, and (2) concurrently it could cause 
nonlinear effects, in addition to (3) the wind effect, described here, making the inner and outer 
arm boundaries overlap because of the different propagating speeds. The dashed lines in Figure [6] 
outlining the maximum density peaks including substructures, a = 2apeak is adopted. 



3.3. Vertical Flatness 

As seen in Figure [2)3, the vertical structure of a gravitational wake has a shape of concentric 
arcs, confined by solid lines representing the empirical formula for the vertical extension limit. The 
angular size of the arcs is best fit to the relation 

^arc = 2tan-i . . , (13) 

which is smaller (or further concentrated t owar d the orbital plane) in comparison to its counterpart 



in a static medium {Mw = 0, dashed; see IK 111 for details). 



In order to quantify the wake mass distribution in the vertical direction, the perturbed density 
a is integrated as a function of latitudinal angle. Figure [7tt displays the integrated wake mass 
normalized by the total wake mass for different wind Mach numbers Myj. The dotted and dashed 
lines (T and SS models) show gradual increases until reaching unity at ^ 50° for J^yj — 3 — 5, 
while the solid lines (SF models) reveal that the saturation is already achieved at ~ 20° for a 
higher wind Mach number Myj ~ 10. Beyond these latitudinal angles, the AGB envelopes remain 
unchanged under the gravitational perturbation of the orbiting substellar object. The decrease of 
the critical angle with implies a more flattened morphology of the circumstellar matter in a 
faster outflow. This is due to the fact that in fast background winds, the secondary's gravitational 
potential relative to the wind kinematics is low so that its perturbation affects only a limited area. 
In slow winds, however, the dependence of the vertical extent on Myj is weak, as revealed in the 
resemblance of the dotted {M.w = 3) and dashed {M.w = 5) lines in Figure EK. 

The distribution of the integrated wake mass is relatively independent of Mp, Cg, Mp, and the 
background wind model. For instance. Figure [7^ shows the overlap of 4 solid lines (SFf, SFm, SFs, 
and SFx models) and 4 dashed lines (SSf, SSm, SSs, and SSx models), respectively, for M.p = 0.5, 
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2.2, 5.0, and 10.0. The independence on the perturber mass Mp and the fluid sound speed Cg (or 
the Bondi accretion radius rs by association) is shown in Figure [7)o (Al, A5-1, and A5-2 models). 
Moreover, the independence on the wind model is found from similarities between the solid lines in 
Figure [7^ (S models) and the lines in Figure [7)3 (A models) for fast winds {A4w — 10) and between 
the dotted (T model) and dashed (S models) lines in Figure [7^ for slower winds (A^iu — 3 — 5). In 
conclusion, the vertical distribution of a wake in an outflowing background is determined by the 
Mach number of the average expansion velocity, (> 1), only weakly depending on the accretion 
and orbital properties of the perturbing object. 



4. DISCUSSION 

In order to probe the existence of planets and/or brown dwarfs around evolved giant stars, 
we investigate the properties of the gravitational wake induced by a substellar mass object in 
an outflowing environment, through hydrodynamic simulations using the FLASH adaptive mesh 
reflnement code. The perturbing object of substellar mass Mp is assumed to be in circular motion 
at distance Vp with the orbital speed of Vp. For an isothermal sound speed Cg, the expansion 
speed Vw{r) of the background wind is determined by the hydrodynamic condition at a reset radius 
representing the wind-generating location. As functions of the object and wind parameters, Mp, 
Vp, Vp, Vw, and Cg, we quantify the observable properties of the wake, revealing that the shape of 
the spiral and arcs are dependent on viewing direction and the density contrast of the structure. 
Conversely, possible future observations of circumstellar spirals and arcs will allow us to place 
constraints on the properties of substellar objects in the winds of AGB stars. 

We find that the overdense arm pattern of the wake propagates outward with the speed of 
V^rm = ^ i Cg, where the upper and lower signs represent the case for the outer and inner 
boundaries of the arm pattern, respectively. Here, the speeds are local, implying that the arm 
spacing varies in the wind acceleration zone and/or in regions where the sound speed varies. In 
order to measure the variation of arm spacing due to the wind accelerat ion, which is believed to 



occur within up to a few tens of stellar radii (jHabing &: OlofssonI |2003| ) . observations with high 



angular resolution (< a few O'.'l for the closest AGB stars) are required. Outside the acceleration 
zone, where the variation of pattern propagation speed is negligible, these pattern propagation 
speeds yield the wake shape in the orbital plane, 

r/rp = ^x{V^± Cs)/Vp + 1 (14) 

for the outer and inner boundaries, respectively, constituting Archimedes spiral shapes. If we define 
the arm spacing Ararm as the interval of the outer or inner boundary individually, the spiral shapes 
of boundaries simply yield 

Ararm = (K« ± Cg) X — ^ (15) 
Vp 

assuming that the pattern propagation speed does not significantly vary within one turn. One can 
estimate the orbital period, 27rrp/Vp, by measuring the observed arm spacing in addition to the 
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wind and sound speeds. Alternatively, this is written by 

_ 2n{V^ ± Cs) 3/2 , . 

from which we can predict the orbital distance rp of the substellar object in the case that the stellar 
mass is known. We note that the sound speed should be taken into account especially in a slow 
wind {Vw ~ Cs). 

The outer arm boundary propagates at a higher speed than the inner boundary by 2cs, leading 
to the broadening of the arm as a function of radial distance. Thus, these high density boundaries 
overlap each other and are in resonance at the overlap positions, rover- Assuming a constant 
wind and sonic speeds, the spacings of these resonant positions is derived to be Arovcr/^^'arm = 
0.5 (Ko/cs T 1) (see eqs. [TT] and [IS])- In cool AGB envelopes, the wind is usually much faster than 
the sound speed iVw ^ Cs), making Arovcr relatively large. For instance, with Vw = lOcg, which is 
easily found for AGB stars, the outer and inner arm boundaries can meet only after five turns. For 
the case of an observational detection of only parts of the spiral, especially when the partial spiral 
is expected over the distance of Arover, a more careful analysis is required to avoid misidentifying 
the outer and inner boundaries. 

We have also provided an empirical formula (eq. [I2]) for the arm-interarm density con- 
trast, a = 5p/pw, along a normalized distance r/vp as a function of V^/cg, V^/cs, and rs/rp = 
GMp/(c^rp). Using this empirical formula, we estimate the properties of a Jupiter wake in the stel- 
lar wind when our Sun bec omes a giant of siz e 1 AU (M* = 0.8 M©, = 2 x 10~'^MQyr~^, 



T^, = SOOO-ftT; according to iHurlev et al.l l2000l ). The wind speed is set to be ~10kms~^ 
b ased on the tren d between the mass-loss rate and the envelope expansion speed (see Fig. 16 
in lFong et al.l[2006l ). and the sonic speed Cg is assumed to be lkms~^. The estimated orbital speed 
Vp = {GM^./rpY^'^ of Jupiter in situ (r^ = 5AU) is 12kms~^, corresponding to the arm-interarm 
density contrast of only 7% of the background density at distance of 100 AU. The density contrast 
does not significantly depend on the orbital distance, i.e., 6-11% for rp = 3-30 AU. Thus, it is 
difficult to detect the gravitational wake of a Jupiter mass object with the current observational 
limitations of sensitivity and angular resolution. In the same range of orbital distance, a 10 Jupiter 
mass object can create a gravitational wake with the density contrast of 30-44 % at 100 AU; and for 
a brown dwarf mass (0.07 Mq), the contrast increases to 160-220 %, which may be detectable with 
the high sensitivity performance of the Atacama Large Millimeter /submillimeter Array (ALMA). 
We note that these numerical values are lower limits since the peak density contrast at the arm 
boundary can be higher with a realistic size of the object much smaller than > 0.1 AU employed 
in this study; but the effect of the object size r<j is not significant unless the size is comparable 
to the accretion radius tb- The required spatial resolution to distinguish the arm pattern separa- 
tion is 12-370 AU depending on the orbital distance [rp = 3-30 AU). On a larger scale, a distance 
corresponding to 5 times this arm separation (in this case of 14,/cs = 10) is the distance between 
overlaps showing a higher density contrast by a factor of 2. 



The case for fast winds compared to the sonic and object speeds is the most relevant case 
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for many observed AGB envelopes (> lOkms"^). For a substellar object orbiting at a distance, 
Tp > 10 {M^/Mq) (14;/10kms~^)~^ AU, the density contrast at large distance (r » Vp, rs) can be 
simply written as 

/rB\i/2 

apeak > y—j for > > Cs, (17) 

from which the measurement of the density contrast provides the value of the Bondi accretion 
radius = GMp/c^ (or the object mass Mp with a given sound speed) for the structure at a 
known distance from the AGB star. Under the fast wind assumption, the density contrast shows 
a shallower decrea se wi th distance, compared to its counterpart in an initially static background 



decreasing as ( Kill ). This shallow decrease in the density contrast facilitates observations at 
great distances from the star, which is favorable to avoid the high obscuration due to either the 
bright stellar glare or the dense central environment. Equation (|17p indicates that the density 
contrast between arm and interarm regions can be greater than Op^ak for distances where 

Assuming that we can distinguish a density fluctuation of 30% {apeak ~ 0.3), the wake due to an 
object of mass Mp would be observable in the central area within a distance of ~ 10 (Mp/Mj) AU 
in the case of Cg = lkms~^. That is, we may expect to observe the wake of a lOMj mass planet 
within 100 AU with a sensitivity corresponding to a density contrast of 30%. Alternatively, it may 
be necessary to take account of the density enhancement due to the planetary wakes when modeling 
the overall profile of circumstellar envelopes using observational data in which the wakes are not 
resolved. 

In the regime for which the background wind is slower than the orbital velocity (Vp ^Vw^Cs), 
corresponding to an object in a close orbit at rp < 10 (M^/Afg) (Ky/lOkms"^)^^ AU, the peak 
density contrast at large distance (r ^ rp, rs) is roughly proportional to the velocity ratio Vw/Vp. 
In this regime, the effect of the wind outflow in opening the spiral pattern competes with the 
object orbital motion tightening the pattern (see eq. [H]), causing a higher density jump at the 
more normal shock boundary to the instantaneous flow. 

Although the analysis in this paper only treats a substellar mass object orbiting about the AGB 
star fixed at the system center, our results can be extrapolated to the gravitational density wake 
of the companion in a stellar binary system. For instance in an AGB envelope, characterized by a 
temperature of 100 K (cs = lkms~^), a solar mass companion would produce a density fluctuation 
of a factor of 1.6 to 1.3 {apeak ~ 0.6-0.3 from eq. [IT]) at distance of 3,000- 10,000 AU from the 



center of the mass losing star, where the spiral pattern is found in AFGL 3068 (jMauron &: Huggins 



20061 ). The density fluctuation from the companion's wake is smaller than that estimated from the 



scattered light observation (the factor of up to ~ 5) , indicating that the spiral pattern of AFGL 306 8 



is likely caused by the reflex motion of the AGB star as sugge sted in iMauron &: Huggind (|2006l ) 



following earlier theoretical work, for example, bv ISokeij (| 19941 ). However, in the central 1" region 



(< 1,000 AU), where the density fluctuation due to the companion is larger, ALMA would possibly 
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resolve a double spiral structure due to the individual stars. To compare the contribution of the 
two different mechanisms, we plan to study the density enhancement due to the reflex motion of the 
mass losing star, similar to the study here. The shapes of two spirals (e.g., pitch angle, arm spacing) 
are expected to be same, because the spiral shape is determined by the local wind properties, Vw 
and Cg. But the spiral of the AGB star is probably detached from the star because of the stellar 
outflow, contrary to the spiral of the companion attached to the object. 

Our work suggests directions for future research. For example, a detailed study of the in- 
ner regions of the AGB envelope is necessary to examine the effects of gas heating/cooling, dust 
anisotropic distribution, and possibly stellar pulsation on the properties of the spiral arm pattern. 
In addition, investigations of binary systems over a wider range of mass ratios should be examined 
to explore the modifications in the spiral/arc patterns in the regime where the effect of the reflex 
motion of the AGB star is also important. Finally, radiative transfer modeling of the molecular line 
emission and dust continuum are encouraged in order to compare theoretical models with observed 
structures in detail. 

We acknowledge a helpful report from an anonymous referee. H.K is grateful to Francisca 
Kemper and Kanak Saha for fruitful comments through reading the manuscript, Paul M. Ricker for 
advice on computational issues, and Noam Soker for discussion on vertical extension. This research 
is supported by the Theoretical Institute for Advanced Research in Astrophysics (TIARA) in the 
Academia Sinica Institute of Astronomy and Astrophysics (ASIAA). The computations presented 
here have been performed through the ASIAA/TIARA computing resource, using FLASKS. code 
developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes 
at the University of Chicago. 
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Table 1. Parameters for numerical simulations 



Object Properties Background Properties Domain Information 



Model 


Mp 






r-p 


Mp 


M, 






Cs 


L 


AL 


refinement 




[Mq] 


[AU] 


[AU] 


[AU] 




[Afoyr-i] 


[AU] 




[kms--"-] 


[AU] 


[AU] 




Accelerated Wind 


via 


Force (MJ*^ = 


OMq) 
















Al 


0.01 


9 


1.0 


10 


5.0 


1 X 10"'' 


2 


11 (10.3) 


1 


300 


0.15 - 2.34 


7-3 


A2 


0.01 


9 


1.0 


10 


5.0 


1 X 10-"* 


2 


11 (10.3) 


1 


300 


0.15 - 2.34 


7-3 


A3 


0.01 


9 


1.0 


10 


5.0 


1 X 10"*^ 


1 


11 (10.3) 


1 


300 


0.15 - 2.34 


7-3 


A4 


0.01 


9 


0.1 


10 


5.0 


1 X 10"® 


2 


11 (10.3) 


1 


300 


0.07 - 2.34 


8-3 


A5 


0.09 


9 


0.1 


10 


5.0 


1 X 10-" 


2 


11 (10.3) 


3 


300 


0.07 - 2.34 


8-3 


A5-1 


0.01 


9 


1.0 


10 


5.0 


1 X 10"" 


2 


11 (10.3) 


3 


300 


0.15 - 2.34 


7-3 


A5-2 


0.09 


9 


1.0 


10 


5.0 


1 X 10"" 


2 


11 (10.3) 


1 


300 


0.15 - 2.34 


7-3 


A6 


0.01 


9 


1.0 


20 


5.0 


1 X 10-6 


2 


11 (10.3) 


1 


600 


0.59 - 4.69 


6-3 


Supersonic Wind Branch (M^^ = 


= 1 Mq) 


















SFf 


0.10 


10 


1.0 


20 


10.0 


1 X 10-s 


10 


10 (9.6) 


3 


1000 


0.98 - 7.81 


6-3 


SFm 


0.10 


10 


1.0 


20 


5J) 


1 X 10-'' 


10 


10 (9.6) 


3 


1000 


0.98 - 7.81 


6-3 


SFs 


0.10 


10 


1.0 


20 


2.2 


1 X 10-'' 


10 


10 (9.6) 


3 


1000 


0.98 - 7.81 


6-3 


SFx 


0.10 


10 


1.0 


20 


05 


1 X 10-" 


10 


10 (9.6) 


3 


1000 


0.98 - 7.81 


6-3 


SSf 


0.10 


10 


1.0 


20 


10.0 


1 X 10-" 


10 


5 (4.2) 


3 


3000 


0.73 - 11.7 


8-4 


SSm 


0.10 


10 


1.0 


20 


5J1 


1 X 10-" 


10 


5 (4.2) 


3 


3000 


0.73 - 11.7 


8-4 


SSs 


0.10 


10 


1.0 


20 


2.2 


1 X 10-" 


10 


5 (4.2) 


3 


3000 


0.73 - 11.7 


8-4 


SSx 


0.10 


10 


1.0 


20 


OS 


1 X 10-" 


10 


5 (4.2) 


3 


3000 


0.73 - 11.7 


8-4 


Transonic Wind Branch (MJ^ = 


1 Mq) 


















T 


0.25 


25 


1.0 


100 


10.0 


5 X IQ--'' 


20 


3.4 (1.7) 


3 


2000 


0.98 - 15.6 


7-3 


Note. 


— Each column represents [1] model 


name, [2] perturbing object's mass Mp 


, [3] Bondi 


accretion 


radius tb = 


GMp/cl [4] 



gravitational softening radius rs, [5] orbital radius rp, [6] orbital Mach number Mp, [7] stellar mass loss rate Mt, [8] stellar size r» that 
generates the wind, [9] wind Mach number Mw on average over the simulation domain (and at rp), [10] sound speed of background gas 
Cs, [11] half of domain size L, [12] range of spatial resolution AL, [13] corresponding range of refinement level with the mother grids of 
64x64x32 in (a:, j/, 2)-directions. Mass and gravitational softening radius of the central star are IMq and 1 AU, respectively. Underlines 
emphasize the main difference of the respective models from the fiducial models appearing at the top in categories "A", "SF", and "SS", 
respectively. 
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Fig. 1. — Background wind Mach number as a function of distance r for (a) transonic hy- 
drodynamic wind (model category "T"), (b) supersonic hydrodynamic wind (model category "S"), 
and (c) wind accelerated by additional force (model category "A"). Thick solid curves show the 
final stable configuration of the background, following the analytic solution (dot-dashed) outside 
the conceptual stellar radius (dotted). Dashed lines in pure hydrodynamic winds indicate the 
critical radius rc = GM^/2cl with constant sound speed Cg, corresponding to the location of the 
sonic radius for the transonic branch in (a) and of minimum velocity for supersonic branch in (b). 
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Fig. 2. — Density enhancement a of model with a fast wind (Al) in logarithmic scale, (a) In 
the orbital plane, the high density region is bounded by quasi-Archimedes spiral arm boundaries, 
r/rp = (/9 X (7W^ lb l)/7Wp + 1 {black dashed and dotted lines) in the polar coordinates (r, (p) with (p 
measured in the clockwise direction. White solid line represents the circular orbit of the perturbing 
object moving along in the counterclockwise direction and currently located at the intersection of 
white dotted lines, (b) In a meridional plane, the wake appears to possess arcs with increasing 
radius along distance from the center. The boundaries are confined by solid lines having angle of 
tan~^[(A1p — 1)^/^(1 + 0.2A4ujA4p)~^] from the orbit. The comparison with the extension limit in 
the absence of wind (i.e., = 0; dashed) demonstrates the tendency for the wake to vertically 
flatten with a fast wind. (A color version of this figure is available in the online journal.) 
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Fig. 3. — Comparison of the shape of spiral pattern between models with the wind Mach number 
Mw of (a) 0, (b) 5, and (c) 10. The perturbing object orbits on a circle at distance {solid, 20 
times greater than the Bondi accretion radius) in the counterclockwise direction with the orbital 
Mach number Mp = 2. Its current location is at the intersection of dotted lines. Dashed line 
denotes the opening angle of the spiral in a static medium, = sin~^A^~^, from which the 
corresponding models with faster winds open up of the spiral pattern. The color bar labels the 
density enhancement a in logarithmic scale. These models are not listed in Table [TJ (A color 
version of this figure is available in the online journal.) 




Fig. 4. — Density enhancement of the spiral structure in logarithmic color scale plotted in the polar 
coordinates (left), in which the slope of the pattern corresponds to the radial propagation speed 
of the spiral arm pattern relative to the orbital speed {right) for T model. Solid and dashed lines 
represent ± 1, tracing the slopes of outer {closed circle) and inner {open circle) boundaries of 
spiral arm, and the vertical dotted lines denote the turns of outer spiral. (A color version of this 
figure is available in the online journal.) 
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Fig. 5. — Density enhancement a for models (a) Al, (b) A2, (c) A3, (d) A4, (e) A5, and (f) A6 
as a function of the distance normahzed by the orbital radius Vp of the perturbing object. The 
colored lines from red to blue display a in the direction apart from the line connecting the object 
and system center with the angle of 0°, 45°, 90°, 135°, 180°, 225°, 270°, and 315° trailing the spiral 
arm. See text for the formula of the complementary black line that outline the peaks. 
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Fig. 6. — Same as Fig. [5] except for models (a) SFs, (b) SFm, (c) SFf, (d) SSs, (e) SSm, and (f) 
SSf. The subpattern showing high density groups of peaks (outHned by dashed hne) is due to the 
overlap between the inner and outer boundaries of the spiral wake pattern. 
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Fig. 7. — Integrated mass M of perturbed density a as a function of the angle from the orbital plane, 
normalized by the total mass enhancement Mtot in the simulation domain, (a) The integrated mass 
of wake increases gradually and levels off at 20° for A4w — 10 cases {solid, SF models), but it reaches 
saturation at 50° for A4w — 3 — 5 cases (dashed, SS models; dotted, T model), (b) Comparison of 
the integrated wake mass in models Al (dashed), A5-1 [solid), and A5-2 [dotted) indicates its weak 
dependence on the accretion radius tb = GMp/c^, and moreover on the perturber mass Mp and 
the sound speed Cg of the gas flow, separately. 



